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We simulate clustering, phase separation and hexatic ordering in a monolayered suspension of 
active squirming disks subject to an attractive Lennard-Jones-like pairwise interaction potential, 
taking hydrodynamic interactions between the particles fully into account. By comparing the hy¬ 
drodynamic case with counterpart simulations for passive and active Brownian particles, we elu¬ 
cidate the relative roles of self-propulsion, interparticle attraction, and hydrodynamic interactions 
in determining clustering and phase behaviour. Even in the presence of an attractive potential, we 
find that hydrodynamic interactions strongly suppress the motility induced phase separation that 
might a priori have been expected in a highly active suspension. Instead, we find only a weak ten¬ 
dency for the particles to form stringlike clusters in this regime. At lower activities we demonstrate 
phase behaviour that is broadly equivalent to that of the counterpart passive system at low tem¬ 
peratures, characterized by regimes of gas-liquid, gas-solid and liquid-solid phase coexistence. 
In this way, we suggest that a dimensionless quantity representing the level of activity relative to 
the strength of attraction plays the role of something like an effective non-equilibrium tempera¬ 
ture, counterpart to the (dimensionless) true thermodynamic temperature in the passive system. 
However there are also some important differences from the equilibrium case, most notably with 
regards the degree of hexatic ordering, which we discuss carefully. 


1 Introduction 

Active matter comprises internal mesoscopic subunits that collec¬ 
tively drive the system far from thermal equilibrium by each in¬ 
dividually converting internal energy into mechanical work. In 
a biological context, examples include crosslinked filaments acti¬ 
vated by molecular motors in the cell cytoskeleton 1 , cells collec¬ 
tively organised in living tissues 2 ! suspensions of motile microor¬ 
ganisms 3 4 , shoals of fish and flocks of birds 5 . Synthetic exam¬ 
ples include vibrated granular monolayers^® and suspensions of 
phoretic colloidal particles® 1 . As has been well documented in 
recent reviews 2K] active systems display a host of exotic emer¬ 
gent effects including swarming, pattern formation, giant number 
fluctuations and nonequilibrium ordering. Among these, this pa¬ 
per focuses in particular on the phenomena of aggregation and 
phase separation in active suspensions. 

In the literature, (at least) two mechanisms have been put for¬ 
ward to suggest that activity should generically lead to aggrega¬ 
tion. The firstEHini operates in systems of active particles that are 
elongated and so have an intrinsic (steric) tendency to align with 
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each other. An activity-mediated coupling between these orien¬ 
tational modes and fluctuations in the local number density then 
provides a mechanism for giant number fluctuations and phase 
separation. The standard deviation AN of particle number in a 
subsystem with a mean number of particles N then scales as N a , 
with a > 1/2. 

The second mechanism 21 does not rely on steric orienta¬ 
tional effects and operates even in systems of spherically sym¬ 
metric particles. It was put forward originally in the context of 
certain species of bacteria that obey simple “run-and-tumble” dy¬ 
namics, in which each particle (at least when in isolation from 
any other particle) swims in straight lines at constant speed be¬ 
tween intermittent tumbling events in which it suddenly changes 
swim direction. A positive feedback loop then arises in which the 
particles (a) accumulate in regions where they move more slowly 
and (b) further slow down in regions where they are impeded 
by crowding from other particles. This renders an initially uni¬ 
form suspension unstable to “motility induced phase separation” 

(Mips)l2U22] 

Following its prediction, MIPS was indeed subse¬ 
quently observed in simulations of spherical (or disklike) active 
particlesISHiZl. 

These simulations however discard hydrodynamic interactions 
between the particles, which are of widespread importance in ac¬ 
tive systems. Motivated by this, in RefJ® W e simulated a suspen- 
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sion of active squirming disks with full hydrodynamics. In do¬ 
ing so we showed that hydrodynamic interactions in fact strongly 
suppress MIPS, with the predicted phase separation replaced by 
only a very weak stringlike clustering. 

We interpreted this suppression in terms of the breakdown of a 
mean-field assumption contained in the (a)-(b) feedback mech¬ 
anism outlined above: if a particle is to experience a slowing 
of its run by crowding from other particles (part (b)), it must 
suffer many collisions during any such run. This assumption 
breaks down in the case of hydrodynamic squirmers, because 
whenever two squirmers scatter off each other hydrodynamically 
they each suffer an 0(1) change in swim direction^. To support 
this interpretation, we further demonstrated suppression of MIPS 
even without hydrodynamics, in active Brownian particles in any 
regime where the time for the decorrelation of swim direction 
fails to comfortably exceed the time between particle collisions. 
The importance of rotational and translational timescales in de¬ 
termining phase separation was also discussed in Refsi^ 31 . 

Nonetheless, strong aggregation effects have been widely ob¬ 
served experimentally in active suspensions such as monolayered 
light-activated colloidal surfers® and phoretic colloidal Janus 
particlesEEH. Conversely, however, many bacterial suspensions 
show pronounced coherent structures in their velocity field with¬ 
out strong spatial inhomogeneities in concentration®®, al¬ 
though aggregation has been studied in suspensions of bacteria 
(with and without active flagellae) subject to an attractive deple¬ 
tion interaction®®. 

In many of these experimental systems, the presence of a slight 
attractive potential between the particles cannot be ruled out. In 
contrast, the simulation study of Refill concerned particles that, 
besides hydrodynamic effects, interact only via a purely repul¬ 
sive potential that onsets steeply at particle-particle contact. The 
original predictions of MIPS likewise concerned purely repulsive 
hard particles 2122 . A key open question, therefore, is whether 
the presence of an attractive contribution to the inter-particle po¬ 
tential might restore a tendency for active suspensions to phase 
separate, even in the presence of hydrodynamics. 

As noted above, the phase behaviour of active particles has 
been extensively studied in the absence of hydrodynamic inter¬ 
actions by means of Brownian dynamics simulations, both with- 
out mf26|30|31|3^ and with El4i 

attractive interactions between 
the particles. In particular, Ref.® reported attraction-dominated 
phase separation at low activity, with a crossover to MIPS at 
higher levels of activity. However these existing studies of at¬ 
tractive active systems again neglect hydrodynamic interactions. 
Given the strong effect of hydrodynamics on the phase behaviour 
of purely repulsive active suspensions 28 , an important outstand¬ 
ing problem is to elucidate the phase behaviour of attractive ac¬ 
tive suspensions with hydrodynamics. 

With that motivation in mind, in this work we simulate a 
suspension of active squirming particles®^ with an attractive 
Lennard-Jones component to the potential of interaction between 
them. Hydrodynamic interactions are taken fully into account us¬ 
ing the simulation method of Ref.®, which allows for short range 
lubrication, long range (power law) propagators, and indeed hy¬ 
drodynamic effects on all lengthscales in between. The method 


moreover gives access to the full range of packing fractions from 
dilute to close packed, thereby allowing a comprehensive study of 
phase behaviour. In these hydrodynamic simulations we neglect 
Brownian motion, considering the limit of purely athermal dy¬ 
namics that is relevant to many active systems, such as bacterial 
suspensions. We return below to justify this neglect more fully by 
commenting on the likely effect of a non-zero temperature. 

To understand the effects of activity in the face of this attractive 
component to the interaction potential, we define a dimensionless 
parameter that quantifies the strength of swimming compared to 
the strength of the attractive interaction. For small values of this 
parameter (low activity), we shall find phase behaviour that is 
dominated by this attractive interaction, as might be expected in¬ 
tuitively. In particular, we find regimes of phase coexistence (gas- 
solid, gas-liquid, liquid-solid, etc.) that closely resemble those of 
passive, thermalised LJ particles. In this sense, the scaled activity 
parameter for the active system might be viewed as a dimension¬ 
less effective temperature that plays a role somewhat analogous 
to the true dimensionless temperature in the counterpart passive 
LJ system. There are however some important differences be¬ 
tween the attraction-dominated phase separation that we report 
here in the active system and its equilibrium counterpart, most 
notably with regards the degree of crystalline ordering. 

In contrast, for high levels of activity we find phase behaviour 
that is relatively unaffected by the LJ attraction. This is also con¬ 
sistent with intuition, because strongly swimming particles are 
expected to be much less affected by the attractive interaction. 
Indeed in this regime we find only a weak tendency for stringlike 
clustering. We give evidence to support the interpretation that 
this is again the signature of a nearby MIPS that has been nar¬ 
rowly avoided by the hydrodynamic suppression mechanism® 
outlined above for purely repulsive particles. 

In between these two limiting regimes - attraction-dominated 
phase separation at low activity, and stringlike clustering at high 
activity - we find a slow crossover region in which the two ef¬ 
fects merge, with a degree of particle clustering that is set by the 
packing fraction and the dimensionless activity parameter. 

Alongside the attractive Lennard Jones contribution to the in¬ 
teraction potential, a steep repulsive interaction is also present 
that diverges at particle contact, i.e., at a particle separation 
r — 2R+, where R is the bare (hydrodynamic) radius of the 
squirmers. We therefore expect much stronger excluded volume 
effects than for the LJ potential that is conventionally used in the 
literature, which diverges only at full overlap r 0 + . In view of 
this, as a warmup to our hydrodynamic simulations we shall first 
map out the phase behaviour of both passive and active Brownian 
particles with this modified LJ potential, in the absence of hydro¬ 
dynamics. This will then provide a benchmark against which the 
results of the hydrodynamic study can be properly compared. 

Because of the heavy computational cost of simulating dense 
suspensions with full hydrodynamics, we have only been able to 
access a relatively small number of particles in each run, typically 
N = 242. Phase diagrams will therefore be elucidated in qualita¬ 
tive outline only, most obviously near any critical point. Indeed 
this provides an additional motivation for first simulating passive 
and active Brownian particles, in order to learn what features 
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of phase behaviour we can realistically expect to capture with 
this small N = 242. Accordingly although we could in principle 
achieve much larger N = 10 5 or even N = 10 6 in the Brownian 
simulations, we use N = 242 here too, to provide a true compari¬ 
son with the hydrodynamic results. While this small N is clearly a 
limitation, it is worth remarking that some of the experiments re¬ 
porting clustering of artificial swimmers also involve a relatively 
small number of particles. We note also that the first simulations 
to map out the phase behaviour of passive Lennard Jonesium used 
a similarly small 

The paper is organized as as follows. In Sec. [ 2 ] we describe the 
Brownian and hydrodynamic models to be studied, along with the 
simulation methods used. Sec. [ 3 ] summarizes the model param¬ 
eters, our choice of units, and the principal adimensional quan¬ 
tities to be explored in the simulations. In Sec. [ 4 ] we discuss the 
initial condition used in the simulations, before in Sec. [ 5 ] defining 
the statistical quantities that we shall measure during the runs. 
In Sec. [6] we present our results for passive Brownian particles, 
active Brownian particles and hydrodynamic squirmers in turn. 
Sec. [7| contains our conclusions. 

2 Models 

In this section we introduce the three models to be studied in the 
remainder of the manuscript. We start by defining the modified 
Lennard-Jones interaction potential that is common to all three 
models. We then define the dynamical rules of the three mod¬ 
els separately: passive Brownian particles (PBP), active Brown¬ 
ian particles (ABP) and hydrodynamic squirmers (HS or simply 
squirmers). Note that each model has its own set of dynamical 
rules, but all are subject to the same modified Lennard-Jones po¬ 
tential. 

For all three models we consider a two dimensional (2D) sys¬ 
tem comprising N disks in a square box of side L with periodic 
boundary conditions. Each particle has a bare radius R (bare 
in a sense to be defined below), implying a bare area fraction 

<j> r = NnR 2 /L 2 . 


2.1 Potential 

The particles experience pairwise Lennard-Jones interactions gov¬ 
erned by the potential 


V L j(r)=4£ LJ 


/OLA 12 /ouy 


V r J 


( 1 ) 


which gives short-ranged repulsion for particle separations r < 
2 1 / 6 olj, and long-ranged attraction for particle separations r > 
2 1 / 6 ct L j. See the dotted line in Fig. [l] Setting <t L j = 2 R then sets 
the minimum of Vlj to be located at r = 2 1 / 6 R. 

As we shall describe in more detail below, the hydrodynami- 
cally squirming disks propel themselves by effecting a tangential 
velocity round their edges, at a hydrodynamic radius R from the 
particle centres. Given the hydrodynamic interactions that we 
shall also describe, short range lubrication effects then in prin¬ 
ciple disallow any particle overlaps in which particles approach 
each other to a separation r < 2 R. However the inevitably finite 
discretisation of time and space in the simulation can in practice 


result in occasional overlaps. To prevent these we follow standard 
practice in the literature and include an additional short ranged 
steep repulsive contribution 


V h (r) = 8 h 


— ( 1 — (r — 2R) — 
r-2Rj \ y ; cr h 


( 2 ) 


for 2 R < r < 2R + <Th, with (r) = 0 otherwise. We then choose 
a shell size = (2 1 / 6 — 1)2 R such that the onset of this steep 
repulsive interaction at particle separations r = 2R + Oh coincides 
with the minimum in the Lennard-Jones potential at r = 2 1 / 6 ct L j = 
2 1 / 6 R. See the thin solid line in Fig. [T] The divergence of this ad¬ 
ditional contribution as r —>> 2R+ then prohibits overlaps in which 
squirmers approach each other to a separation of less than twice 
their hydrodynamic radius. 

For simplicity we set the potential amplitudes equal through¬ 
out, £ LJ = e h = £, and use a value for the exponent v = 2. The final 
combined pairwise interaction potential V = V LJ 4 - is shown by 
the thicker solid line in Fig. [l] Although the additional repulsive 
shell provided by Vh is not needed for the PBP or ABP (recall that 
we include it simply to avoid occasional hydrodynamic overlaps 
in the squirmers), we use it for the Brownian models too to en¬ 
sure a consistent potential in all simulations, thereby allowing a 
truly direct comparison between all three models. 

The presence of this steeply repulsive shell acts to confer an in¬ 
creased effective disk radius R e ff = R + o\/2 that slightly exceeds 
the bare disk radius R. Accordingly we define an effective area 
fraction 0 = NnR 2 ^/L 2 that likewise slightly exceeds the bare 
one. Whenever we quote an area fraction throughout the rest of 
the paper, it is this effective quantity 0. 
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Fig. 1 Modified Lennard-Jones potential in units of e=1, R= 1. 
Lengthscales cr L j = 2 and cr h = 2 ( 2 h 6 - 1 ). 


2.2 Dynamics 

2.2.1 Passive Brownian particles (PBP). 

Before considering active suspensions, which are the primary fo¬ 
cus of this work, we will first consider passive Brownian particles 
(PBP) subject to the same modified Lennard-Jones potential dis¬ 
cussed above. We shall do so in order first to establish the equilib¬ 
rium phase behaviour of a system of particles with this modified 


This journal is ©The Royal Society of Chemistry [year] 


Journal Name, [year], [vol.], l|||3 












potential, which will then serve as a point of reference for study¬ 
ing the non-equilibrium phase behaviour of active systems subject 
to the same potential, both without (ABP) and with (HS) hydro- 
dynamic interactions. 

Within the model of passive Brownian dynamics, the position 
r i of the zth particle evolves according to 


dxi 

dt 


]_ dV(|r/ r j|) 

St, 


+ V2DZi, 


(3) 


tides. The equivalent rotational diffusion coefficient that char¬ 
acterises these active phenomena is expected be much larger in 
magnitude than the true thermal one. In order to focus on the ef¬ 
fects of activity, therefore, we have removed the thermal noise 
from the translational equation altogether, setting the transla¬ 
tional diffusion coefficient D = 0 upfront® 25 . Indeed we suggest 
that use of the terminology “Brownian” for ABPs is perhaps bet¬ 
ter avoided altogether, though we keep it for consistency with the 
earlier literature. 


with the effect of a solvent coarse-grained in terms of an effective 
drag of coefficient y, and with a thermal noise characterised by a 
diffusion coefficient D. The Gaussian random variable S has 


<s i,o(()> = 0, 


<S(0i,aEy,p(O> = Sij5 a p5{t-t'), (4) 


where Latin and Greek symbols denote particle number and 
Cartesian coordinate respectively. The drag is related to the dif¬ 
fusion constant by the Einstein relation D = k^T /y, where T is 
the thermodynamic temperature and &b is Boltzmann’s constant. 
The particles also suffer stochastic rotation following Eqn.|6] be¬ 
low, with D v = 3D/4R 2 . However this is irrelevant, because the 
translational dynamics is independent of particle orientation or 
rotation for these passive particles. 


2.2.2 Active Brownian particles (ABP). 


Within the model of active Brownian particles, the translational 
velocity of the zth particle is given by 


dxi 

dt 


1 y dV(|r,--r;|) 
dr ‘ 


+ v fl e/. 


(5) 


Besides the effect of the conservative force stemming from the 
modified Lennard-Jones potential, each particle also moves with a 
self-propelling velocity v a e*. The speed v a is an imposed constant 
of the problem. The unit vector e* = (cosa*, sina*-) prescribes the 
instantaneous preferred direction of self-propulsion, and obeys 
stochastic angular dynamics governed by a rotational diffusion 
constant D r , such that 

~ = v'3A'S,■(./). (6) 

The Gaussian random variable E, obeys 


< S i(t) > = 0, 

<S i(t)2j(t')> = 8ij8(t-t'). (7) 

Note that in an equilibrium system the rotational diffusion co¬ 
efficient D v would have to relate to the corresponding transla¬ 
tional diffusion coefficient D as D r = 3D/4R 2 . However the ac¬ 
tive Brownian particles of interest here do not need to follow any 
such relation. Indeed, the stochastic angular dynamics prescribed 
above is not intended to be of true thermal original, but rather 
as a continuous-time statistical model of (for example) tumbling 
events in a system of run-and-tumble bacteria, or of particle re¬ 
orientations occurring via hydrodynamic scattering off other par- 


2.2.3 Hydrodynamic squirmers (HS). 

The model of hydrodynamic squirmers is based on a minimal de¬ 
scription of microbial propulsion intended to mimic locomotion 
through a fluid by the beating of many cilia on the surface of a 
microbe | 49 | 50|54 | 55 |64| jYie particles are assumed to be neutrally 
buoyant and so force free. Their size and swimming speed are 
sufficiently small that the Reynolds number is assumed negligi¬ 
ble (zero inertia) and their swimming dynamics purely Stokesian. 
Their size is however large enough that Brownian motion (both 
rotational and translational) is assumed negligible, correspond¬ 
ing to an infinite Peclet number. Other sources of rotational noise 
such bacteria tumbling are also neglected, and we return to justify 
this assumption carefully below. 

Consider then a monodisperse suspension of N inertialess disk¬ 
like swimmers, each of radius R, and each propelling itself 
through an incompressible Newtonian fluid in the x — y plane by 
means of an active “squirming” motion modeled by a prescribed 
tangential velocity 

S(6 — apj = B\ sin(0 — af) + sin2(0 — a*) (8) 

round the disk edge, at a distance R from the disk’s centre. For 
simplicity this tangential velocity is taken to be time-independent, 
representing an average over many beating cycles. We denote by 
P = B 2 /B 1 the ratio of the second and first modes. The sign of /3 
then distinguishes between “pushers” and “pullers”, with positive 
P corresponding to pullers and negative P to pushers. The am¬ 
plitude of P prescribes the degree of polarity, with |/3| <C 1 giving 
strongly polar swimmers and |/3| > 1 apolar “shakers”. 

A single disk undisturbed by any others in an infinite box 
then has a swim speed vo = B i/2 54 , with an instantaneous pre¬ 
ferred swim direction e* = (cos a*, sin a*). In a suspension of many 
disks the actual swim speeds and directions evolve over time due 
to hydrodynamic interactions mediated by the Newtonian fluid 
through which the particles move, as we now describe. 

The velocity and pressure fields v(r,f) and p(r,t) in the fluid 
outside the disks obey the condition of mass balance for incom¬ 
pressible flow 

0 = V-v(r,0, (9) 

and the condition of force balance in the inertialess limit 

0 = i]V 2 v(r,z) - V/?(r,f) +f. (10) 

Here rj is the fluid viscosity and r = (jc,y) is a position vector in the 
v —y plane. (For convenience the interior of the discs is also taken 
to contain Newtonian fluid obeying Eqns. [ 9 ] and [lO] We then 


4 | Journal Name, [year], [vol.],i H 


This journal is ©The Royal Society of Chemistry [year] 








simply discard this part of the mathematical solution, because we 
are only interested in the behaviour of the fluid outside the disks.) 

We recognise Eqn. [To] as the Stokes equation subject to addi¬ 
tional source terms f. These represent forces round the edge of 
each disk: 

f(r,;)=l>(e«)S(r,--fl), (11) 

i 

in which we are summing here over the separate polar coordinate 
systems (r*, Of) of the disks, such that for the /th disk 

r = R/(f) + r/cos(0;)x + r; sin(0/)y, (12) 

where R/(t) is the position of the disk’s centre. 

These forces are included so as to ensure that the /th disk has 
at any instant a velocity round its edge 

Vf = Vi-R£ii6i + S(e - ai)Gi , (13) 

comprising solid body translation and rotation, plus the tangen¬ 
tial squirming motion prescribed by the slip velocity function in 
Eqn. [8] above. 

At any timestep t—>t+Dti n the simulation, given a current set 
of particle positions R; and preferred swim directions known 
quantities are the lowest two modes of the force (which are con¬ 
strained to ensure zero total force and torque for each disc, con¬ 
sistent with these particles being swimmers driven by their own 
internal dynamics and not subject to externally imposed force 
monopoles or torques), as well as all higher modes of the veloc¬ 
ity (as prescribed by the tangential swimming function 5). Un¬ 
knowns to be calculated are the lowest two modes of the ve¬ 
locity, giving the solid body translational velocity V* and rota¬ 
tional speed (Also calculated as part of this process are all 
higher modes of the force, but we do not report these). The par¬ 
ticle positions and swim directions are then updated according to 
R ( —y R i T Dt\i and ccf — y oci T Dt£li. 

In this way we exactly (once the numerical mode numbers have 
been converged upon) solve for the Stokes solution in the fluid 
outside all the disks, subject to the boundary conditions of the 
prescribed tangential velocities on the disk edges. Out of this cal¬ 
culation emerge all effects of hydrodynamic interaction, includ¬ 
ing long range power law propagators, short ranged lubrication 
effects, and the physics on all lengthscales in between. 

For active particles, the Peclet number Pe is defined as the time 
taken for a particle thermally to diffuse a distance equal to its own 
radius, divided by the time taken for it to swim the same distance. 
Here we have assumed this to be infinite, as in Refs.^3 23 sup¬ 
pressing Brownian motion entirely and considering only the de¬ 
terministic hydrodynamics defined above. This is the relevant 
limit physically for many active suspensions of, e.g., swimming 
bacteria. We do, however, return to comment further on Peclet 
number below. 

With less justification a priori, we also ignore any intrinsic 
tumbling dynamics of the particles, such that within our model 
angular reorientation occurs only via hydrodynamic interactions 
as the particles scatter off each other. However, as we will ex¬ 
plain further in the results section below, we expect this neglect 
of tumbling to be unimportant, because at the high volume frac¬ 


tions of interest here we find that the relevant physics comes from 
the high reorientation rate set by these scattering events, which 
would only be increased still further by the explicit inclusion of 
tumbling events. 

In the literature both spherical (3D)E3 and cylindrical ( 2 D)EH 
versions of this squirmer model have been studied. For computa¬ 
tional efficiency we use the 2D model®, which in fluid mechani¬ 
cal terms effectively corresponds to infinitely long (in the z direc¬ 
tion) squirming cylinders propelling themselves in the plane v — y 
that is orthogonal to their length. While this geometry is of course 
to some extent artificial (as in all studies of reduced dimension¬ 
ality), it is important to realise that concerns about the famous 
Stokes catastrophe of 2D hydrodynamics are not relevant here. 
That catastrophe applies specifically to the case of a particle sub¬ 
ject to a force monopole, whereas the swimmers of interest here 
are force free overall, experiencing only higher order force multi¬ 
poles. Indeed in a previous publication^ we carefully compared 
this 2D case of infinite cylinders to a more realistic 3D scenario of 
squirming disks of finite thickness (i.e. highly flattened cylinders) 
in a highly viscous film of the same thickness, surrounded by a 
fluid of much lower viscosity. Reassuringly we found no quali¬ 
tative difference between these cases, giving confidence that our 
2D calculations provide a good representation of the physics of 
suspensions of squirmers with hydrodynamics. 

3 Units and parameter values 

Table [l] summarizes the parameters of the problem, defining the 
simulation geometry, particle size, modified Lennard-Jones poten¬ 
tial and the dynamics of the three models. Throughout the paper 
we work in units of length in which the particle radius R = 1, 
of modulus in which the potential energy £ — 1, and of time in 
which the drag coefficient y = 1 (for Brownian particles) or the 
Newtonian viscosity q = 1 (for hydrodynamic squirmers). 

Because of the heavy computational cost of the hydrodynamic 
simulations (recall that we solve the Stokes equation fully in the 
plane, recovering out of this calculation both far field power law 
propagators, near field lubrication effects, and the physics on all 
lengthscales between these), the number of squirmers that we 
can feasibly simulate is strongly constrained. Here we choose 
N = 242. The more obvious choice of N = 256 carries a greater 
risk of lattice defects at area fractions approaching close packing. 
We have also checked that the phenomena we report are robust 
against doubling N. For the Brownian dynamics (both passive 
and active) many more particles would be possible because of the 
much lower computational demand of those simulations: up to 
N = 10 5 or N = 10 6 Brownian particles can be simulated in run 
times comparable to those for N = 242 hydrodynamic squirmers. 
However we restrict ourselves to N = 242 for the PBP and ABP as 
well, in order to ensure a direct comparison between all the mod¬ 
els. We return at the start of Sec. [6] below to comment carefully 
on the extent to which one might expect to capture a system’s 
phase behaviour with such a small number of particles. 

The lengthscales Olj and <Jh associated with the potential are 
set as discussed in Sec. |2,1| above and summarised again in ta¬ 
ble □ 

An important parameter that then remains to be varied in the 
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Quantity 

Symbol 

Dimensions 

Dimensionless form 

Value 

PBP 

ABP 

HS 

Particle radius 

R 

L 

1 

length unit 

/ 

/ 

/ 

Potential energy 

£ 

GL 3 

1 

modulus unit 

/ 

/ 

/ 

Drag coefficient 

7 

GTL 

1 

time unit for PBP/ABP 

/ 

/ 

X 

Solvent viscosity 

n 

GT 

1 

time unit for HS 

X 

X 

/ 

Number of particles 

N 

1 

N 

242 

/ 

/ 

/ 

Potential shell size 

Oh 

L 

dh = f 

2 (2*/ 6 — 1) 

/ 

/ 

/ 

Lennard-Jones length 

°lj 

L 

d U=*¥ 

2 

/ 

/ 

/ 

Effective area fraction 

0 

1 

0 

to be varied (all models) 

/ 

/ 

/ 

Thermal diffusion coefficient 

D 

L 2 T~ 1 

D=^ 
r-\ _ AA 

b V a 

to be varied (PBP) 

/ 

X 

X 

Rotational diffusion coefficient 

D v 

T~ 1 

to be varied (ABP) 

- 

/ 

X 

Active swim speed 

V a 

LT~ l 

Va 

to be varied (ABP + HS) 

X 

v a yR 

£ 

B\r\R 2 

2e 

Active polarity 

P 

1 

P 

to be varied (HS) 

X 

X 

/ 


Table 1 Parameters for passive Brownian particles (PBP), active Brownian particles (ABP) and hydrodynamic squirmers (HS) with modified 
Lennard-Jones interactions. Note that PBP are rotationally symmetric so rotational diffusion is unimportant. We work in dimensions space GLT of 
modulus G, length L and time T, rather than the more usual MLT with M being mass. 


simulations for all three models is the area fraction 0. Because of 
the strong repulsive force characterized by the small shell length 
Oh, it is reasonable to work with the effective area fraction defined 
above rather than the bare area fraction 0 r = NnR 2 /L 2 . The real 
area fraction in practice of course depend on the ratio between 
the temperature (or active velocity) and the potential energy e. 
However the steepness of the repulsive shell’s potential renders 
this effect small in practice. Therefore in what follows, unless 
stated explicitly otherwise, we report the effective area fraction. 
Our simulation method allows us to study high area fractions even 
for hydrodynamic squirmers, and accordingly we shall explore 
right across the range from dilute towards close packed. 

Besides the area fraction, the other important parameters that 
remain to be varied are as summarised follows. For the passive 
Brownian particles we have the adimensionalised translational 
diffusional coefficient D , which is also equivalent to the adimen¬ 
sionalised thermodynamic temperature. For the active Brownian 
particles we have the adimensionalised active swim speed v a and 
adimensionalised rotational diffusion coefficient £ _1 . Finally for 
the hydrodynamic squirmers we have the adimensionalised active 
swim speed v a and the stresslet parameter /3. We shall return to 
discuss in more detail the physical meaning of each of these adi- 
mensional quantities in each of the individual subsections con¬ 
cerning PBP, ABP and HS below. 

4 Initial conditions 

In each run the system is initialised by placing the particles at ran¬ 
dom positions and with orientations likewise chosen at random 
from a top hat distribution between 0 and In. Obviously some 
particles will inevitably overlap after this procedure. Overlaps 
are then removed by evolving the system with zero-temperature 
dynamics in a landscape of mutual particle repulsion. 

We anticipate (but have not checked explicitly) that in any re¬ 
gions of the phase diagram where our simulations attain steady 
state, that state would be independent of the initial conditions. 
However, as we shall discuss below, at low temperatures and area 
fractions we typically find a regime of coarsening dynamics in 
which the average cluster size slowly increases with time. Due 
to the very slow dynamics we have not been able to attain steady 
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state in any feasible run time in that regime, and these states must 
inevitably depend on the initial conditions. 

5 Measured quantities 

In this section we describe the various quantities that we mea¬ 
sure in our simulations in order to analyse clustering and phase 
behaviour across the different models. 

5.1 Snapshots 

In each case we first make a simple visual observation of snap¬ 
shots of the system’s configuration, on a grid of values of the area 
fraction and adimensionalised temperature (or activity parame¬ 
ter). Unless otherwise stated the snapshots shown are in steady 
state, verified by confirming that the time signals of the mean 
cluster size, mean particle velocity, etc. have attained their fi¬ 
nal values (which are of course subject to fluctuations). Where 
stated, in some cases of extremely slow coarsening at low tem¬ 
perature, the snapshots shown are for the longest times that are 
feasibly accessible numerically. 

5.2 Cluster analysis 

We statistically analyse the degree of clustering in the system us¬ 
ing a single-link algorithm, the details of which are described in 
Refs.^3 57 . Any two particles are taken to be in the same cluster 
if their separation r is smaller than a cutoff value r = r c . Indeed it 
is worth reiterating that any two groups of particles are taken to 
belong to the same cluster even if only a single pair of particles, 
one in each group, satisfies this cutoff condition. We choose a 
value for the cutoff separation r c = 2.4 R, which just exceeds the 
value of the particle separation at the minimum in the potential 
at r = 2 l ^2R w 2.24 R. We have checked the robustness of our re¬ 
sults to reasonable variations in the value of this cutoff. Common 
practice in the literature ^—I is to use a cutoff smaller than the 
first minimum of the radial distribution function, and we have 
followed that convention here. 

Using this algorithm we measure the cluster size distribution 
p(n,t) at any instant of time t, and the first moment h(t) of 
this distribution, which gives the mean number of particles per 
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cluster. We usually report this quantity once the system has 
reached a statistically steady state after many time units, and 
time-averaged over a time interval At to ensure good statistics. 
Typically At = 310 3 R 2 /3D for PBP, At = 3 x \0 3 R/3v a for ABP and 
At = !0 3 R/3v a for HS. In a few, explicitly stated, cases of extremely 
slow coarsening at low temperature the cluster size is reported at 
the longest times that are feasibly accessible numerically, again 
time averaged over At. 

The instantaneous (non time averaged) behaviour of n(t) as 
function of time t (after initialising the system in a disordered 
configuration at time t = 0 as described in Sec. B will also be 
used to to characterise the dynamical coarsening of clusters in 
the system. 

5.3 Particle number fluctuations 

In order to gain insight into fluctuations in the local area fraction 
of particles, we divide the simulation box into a grid of G x G cells 
each of dimension Lq = L/G, with a typical value of G = 5. We 
then define the instantaneous local particle number in each cell i 
of the grid as Ni(t) and write this as the sum of the average value 
plus a fluctuation away from this: 

Ni(t)=^+8Ni(t), (14) 

and measure the number fluctuations as 

§ N = (y/(Shh(t) 2 )i)„ (15) 

which we typically report scaled by VN, with N= (Ni(t))i = N/G 2 . 
Here (-) t represents an average over time, and (•),• an average over 
particles. 

Alternatively we can think in terms of the local area fraction 
(f>i(t) =Ni{t)nR 2 /Lq = Ni(t)(j)G 2 /N, which for a given grid division 
G and total number of particles N is related to the local number of 
particles N((t) by a constant factor that depends only on the global 
average area fraction 0. The local area fraction fluctuations are 
then defined as 

8i = {y/{8m 2 )i)t- (16) 

To connect these two measures we note simply that 8^/8 = 
G 2 8n/N. In what follows below we sometimes find it convenient 
to report <5/v, which more effectively shows fluctuations at low 0, 
and sometimes 8^, when focusing on higher 0. (The values of 
G = 5 and N = 242 are the same in all our simulations.) 

5.4 Hexatic order parameter 

While the cluster size distribution and particle number fluctuation 
measures described so far are effective at discerning the degree 
of aggregation, they are poorly suited to differentiating between 
liquid and solid phases. We therefore define a hexatic order pa¬ 
rameter T'g to do this. 

We start by defining an instantaneous particle-based local hex¬ 
atic order, for a given particle j at any time t: 

f'6;(0 = 4 (17) 

n j k= 1 


Here hj is the number of neighbours of particle j, defined as those 
particles inside a threshold separation r <r c = 2.4 R, while Qjk(t) 
is the angle of the bond between particles j and k. 

From this instantaneous particle-based measure of local hexatic 
order one can then define various different measures of global 
hexatic order. Seja for two definitions. Here we focus specifi¬ 
cally on the one that averages in time (-) t first for each parti¬ 
cle, and then averages the modulus of this quantity over all parti¬ 
cles (-)j to give 

'Pfi = dCMOM);- CIS) 

For a single perfectly hexatically ordered crystal we expect = 1, 
while for a totally disordered phase T'g = 0. 

6 Results 

We now present the results of our simulations. Even though the 
equilibrium phase behaviour of passive Lennard-Jones (LJ) par¬ 
ticles has been mapped out previously, both in 2 dE 3 and 3 dE 3 
we shall nonetheless start in Sec. |6.1| by studying the equilibrium 
phase behaviour with our modified LJ potential, which includes 
the repulsive shell added to avoid particle overlaps for the squirm¬ 
ers. We do so in order to allow in later sections [672] and [673] a truly 
direct comparison of the non-equilibrium phase behaviour of ac¬ 
tive particles - both without and with hydrodynamics - with the 
equilibrium phase behaviour of their passive counterpart experi¬ 
encing exactly the same interaction potential, including this re¬ 
pulsive shell. (We do not perform simulations of passive particles 
with hydrodynamics, because the phase behaviour of equilibrium 
systems depends only on the static interaction potential and not 
the dynamics.) 

Because the number of particles that we can access for the hy¬ 
drodynamic squirmers is heavily limited by computational cost 
(we use N = 242), it is furthermore important to establish in this 
more familiar equilibrium context to what extent phase behaviour 
can feasibly be discerned with such a small system size. Although 
this number of particles is of course very small for systems with¬ 
out hydrodynamics, by the standards of today’s computational ca¬ 
pabilities, it is worth noting that the early papers to map out the 
phase behaviour of Lennard-Jones and hard sphere systems did so 
using typically 200—1000 particles 51 53 , even in 3D. Indeed these 
early works established that surprisingly small numbers of parti¬ 
cles are needed to map the overall features of the phase diagram 
correctly. Accordingly we claim that our N = 242 simulations can 
capture the main features of phase behaviour. However we do not 
claim full quantitative statistical analysis of two phase regions, or 
even less of critical points, to be a realistic goal in this work. 

6.1 Passive Brownian particles 

As an aid to the reader’s memory we start by sketching in Fig. [2] 
the main features of the equilibrium phase diagram of particles 
experiencing unmodified LJ interactions, in the plane of adimen- 
sionalisation diffusion coefficient D versus the area (or volume) 
fraction 0. Recall that D = Dy/s = k^T/e measures the thermal 
energy of the system relative to the characteristic energy of the 
modified Lennard-Jones potential. As noted above this phase dia- 
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Fig. 3 Passive Brownian particles. Snapshots of the system’s configuration on a grid of values of the scaled thermal diffusion coefficient D and area 
fraction </>. Each snapshot is taken at a long time t = 10 3 R 2 /D after the system was initialised at t = 0 in a random state as described in Sec. [ 4 ] Note 
that the scale is nonlinear at the largest D and at low 0. 
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Fig. 2 Sketch of the main qualitative features of the equilibrium phase 
diagram of disks with a Lennard-Jones interaction potential VuEHm. 
S=solid, L=liquid, G=gas. 

gram has been comprehensively mapped out previously 51 52 . In¬ 
dicated are the single phase regions of gas (G), liquid (L) and 
solid (S), and the phase coexistence regions G-L, G-S, L-S. We 
emphasize that Fig. [ 2 ] is a rough qualitative sketch only with no 
attempt at quantitative accuracy 

Turning now to our modified Lennard-Jones potential, we show 
in Fig. [3] snapshots of the system’s configuration on a grid of 
values D and area fraction 0. All snapshots shown are for the 
longest time that is feasibly accessible computationally. At low 
temperatures ,D ~ 0.4, aggregated domains are evident across the 
full range of area fractions. These coarsen very slowly as a func¬ 
tion of time, and at the lowest area fractions and temperatures 
(typically in the rectangle of phase space with 0 < 0.125 and 
D < 0.4) this process cannot be evolved fully to steady state within 
a feasible run time. (For example, an order of magnitude esti¬ 
mate suggests that for the simulation at 0 = 0.0125 and D = 0.05 
to have fully coarsened, a run time of about 6 months involving 
O(10 9 ) timesteps would have been needed.) The snapshots in this 
region are therefore not in true steady state, though we are confi¬ 
dent that the qualitative features will not change even for longer 
run times, which would simply result in slightly larger clusters. At 
higher temperature we have checked that the system is in steady 
state. 

Simple visual inspection of these snapshots suggests an overall 
structure to the phase diagram that indeed captures the main fea¬ 
tures of the sketch in Fig. [ 2 ] In particular, a G-S coexistence region 
is visually apparent for D less then about 0.35. Above this, a G-L 
coexistence region extends up towards what resembles a critical 
point somewhere in the region of 0 = 0.4 and D = 0.55. However 
we do not focus on this critical point in detail because we do not 
have sufficient particles in our simulation to characterise it prop¬ 
erly. For high area fractions, roughly 0 > 0.7, solid-like ordering 
is evident even at high temperatures. This area fraction for the 
onset of solid ordering is lower than that expected for particles 
with conventional LJ interactions, but is consistent with the fact 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

0 

(a) 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 


0 

(b) 



0 

(c) 


Fig. 4 Phase diagram of passive Brownian particles mapped using (a) 
normalised mean number of particles per cluster n* = n/N, (b) particle 
number fluctuations 8 n /Vn for G = 5 and c) hexatic order parameter T^. 


that our modified LJ potential includes an effectively hard shell 
repulsion at r = 2 R. 

To try map the phase diagram more quantitatively, we show in 
Fig.|4]colour maps of the statistical measures introduced in Sec. [5] 
above. Subfigure (a) shows the mean number of particles per 
cluster, normalized by the total number of particles in the system: 
ft = h/N. Note that a system-spanning homogeneous phase with 
all particles in same cluster would have ft = 1, while a gas phase 
with no clusters would have ft = 1 /TV, which is small for large N. 
With this in mind we have also marked on the colourmap isolines 
of ft = 1/50,1/10,1/3, which correspond to n = 4.84,24.2,80.7 
for our N = 242. We suggest that the region in between the first 
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Fig. 5 Representative isolines of particle number fluctuations 
8 n /Vn = 1.5 and ¥ 6 = 0.1,0.3. Regions broadly consistent with the 
sketch of Fig[2]are evident. 

two of these corresponds to one of clusters of a finite size n> 1, 
but which do not yet span the full system size n* = 1. Compar¬ 
ison with simulations of larger N would however be needed to 
substantiate this claim properly. Note that, as above for the snap¬ 
shots, in the region roughly given by 0 < 0.25 and D < 0.4 the 
clusters are still coarsening so that waiting for an even longer 
run-time would result in a larger ft. 

Subfigure (b) shows the particle number fluctuations 8^/VN, 
measured by dividing the simulation box into a grid of G x G = 
5x5 subunits. Clearly apparent is a region with strongly en¬ 
hanced number fluctuations, which we suggest can be interpreted 
as a region of phase coexistence. However this measure provides 
no means of distinguishing between a phase coexistence that in¬ 
volves only disordered phases, G — L, from one involving the or¬ 
dered solid phase, G - S or L - S. In any case it is strongly sup¬ 
pressed at high volume fractions 0 > 0.7 and so fails altogether as 
an indicator of phase separation in this regime. 

To try identify these solid regimes, we show in subfigure (c) a 
colour map of the hexatic order parameter. The bright region for 
0 > 0.75 clearly indicates the presence of crystalline ordering. We 
further suggest the region D < 0.35 to be one of G — S coexistence 
(or fully solid at high 0, but this cannot be discerned separately), 
consistent with visual inspection of the snapshots in Fig. [3] (For 
small area fractions the solid component is too weak to show in 
this representation. As noted above, however, in this part of the 
phase diagram coarsening is so slow that this system does not 
attain a true steady state within the run time.) We also suggest 
that for D > 0.35 the region 0.75 < 0 < 0.8 is one of L — S coex¬ 
istence, and that of 0 > 0.8 a homogeneous S phase. We empha¬ 
size however that these boundaries are suggestive only, and that 
simulations with far more particles would be needed to provide 
quantitatively conclusive evidence. 

Finally we assemble in Fig. [5] a representative isoline of the 
particle number fluctuations 8^/Vn =1.5 together with two rep¬ 
resentative isolines of the hexatic order parameter T'g = 0.1,0.3, 
revealing a phase diagram that indeed qualitatively resembles the 
sketch in Fig. [2] 

Although our main aim in this section has been to report equi¬ 


librium phase behaviour, we return finally to reiterate that in the 
rectangle of phase space with 0 < 0.25 and D < 0.4 the system 
doesn’t equilibrate in any feasible run time but rather shows a 
slow coarsening behaviour in which the typical cluster size in¬ 
creases as a power law n ~ ft. A value for the exponent /3 = 2/3 
provides an excellent fit to the data (not shown), consistent with 
the expectation® of a linear domain size growing as L ~ r 1 / 3 . Al¬ 
though we have been careful only to extract data for the coarsen¬ 
ing exponent in regimes where the typical cluster size is less than 
the system size, a more detailed study of the model’s coarsen¬ 
ing dynamics, for which finite size effects should be investigated 
carefully, is out of the scope of this paper. 

6.2 Active Brownian Particles 

In the previous subsection we mapped out the equilibrium phase 
behaviour of passive Brownian particles (PBPs) in the plane of 
area fraction 0 and adimensional diffusion coefficient D = Dy/s = 
k^T/s, which measures the system’s thermal energy compared 
with the characteristic energy e of the modified Lennard-Jones 
potential. As discussed, we recovered the main features of phase 
behaviour expected from previous studies of Lennard-Jones par¬ 
ticles, with some differences consistent with our use of an addi¬ 
tional steep repulsive contribution Vh- 

We turn now to active Brownian particles (ABPs). Although 
the phase behaviour of ABPs with unmodified Lennard-Jones in¬ 
teractions has been studied previously®, we will map it out here 
for our modified Lennard-Jones potential using the same small 
number of particles N = 242 as will be feasible to explore compu¬ 
tationally for the hydrodynamic squirmers in Sec. |6.3| below. This 
will then enable a truly direct comparison between ABPs, which 
lack hydrodynamic interactions, and squirmers, which have them. 

As discussed above, the ABPs move via ballistic swimming at 
speed v a combined with stochastic reorientational dynamics with 
angular diffusion coefficient D v . The latter is intended as a con¬ 
tinuous time model of (for example) discrete tumbling events, in 
some species of bacteria, or reorientations due to hydrodynamic 
scattering off other swimmers. It is not intended to model true 
thermal angular diffusion, which would be much smaller in com¬ 
parison. Indeed, as noted above, we consider here the athermal 
limit in which any thermal contribution to D v is zero, which is the 
relevant limit for many swimming microbes. 

Despite this lack of true thermal translational diffusion, an ob¬ 
vious question to ask is whether the activity of the ABPs (ballistic 
swimming combined with stochastic reorientation) plays a role 
somewhat akin to the thermal diffusive dynamics of the PBPs, 
such that the non-equilibrium phase behaviour of the ABPs is in 
some degree analogous to the equilibrium phase behaviour of the 
PBPs. For example, we might intuitively anticipate phase separa¬ 
tion at small v a due to the attractive effects of the Lennard-Jones 
potential, with a return to homogeneous behaviour at large v a 
where the activity is strong enough to overcome that attraction. 
In what follows, therefore, we will seek to compare the equilib¬ 
rium phase behaviour of PBPs in the plane of (5,0), as mapped 
out above, with the non-equilibrium phase behaviour of the ABPs 
in the plane of (v a , 0), with v a first suitably adimensionalised. 
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However, whereas for PBPs the only two relevant parameters 
were indeed D and 0, for ABPs the rotational diffusion coeffi¬ 
cient D x is also important, alongside v a and 0. The reason for this 
is as follows. In the case of spherical PBPs rotational diffusion, 
though surely present physically, is irrelevant: each particle is 
rotationally symmetric and the translational dynamics is accord¬ 
ingly independent of particle orientation. In contrast, for ABPs the 
translational dynamics depends strongly on the orientational dy¬ 
namics, because the instantaneous preferred swimming direction 
of any particle is (by definition) prescribed by its instantaneous 
orientation. In view of this we must actually consider the phase 
behaviour of ABPs in the plane of (v a , 0) for several different val¬ 
ues of D r , with both v a and D r first suitably adimensionalised. 

To adimensionalise the rotational diffusion coefficient D v we 

define a parameter®! 


r = Va ?o 

^ _ D r Rf(<j>) ~ T C (0)’ 


(19) 


which is the ratio of the characteristic decorrelation time of the 
particle’s orientation (swim direction) 


t 0 = 


1 


A 


( 20 ) 


with the characteristic time interval between particle collisions 


T c = 


RfW 

V a 


( 21 ) 


Here /(0) is an adimensional function that we expect to de¬ 
crease with increasing 0. In the dilute gas limit of small 0 we 
would estimate / = tt/0 from mean free path arguments, while 
for more crowded systems with 0 = 0(1) we might instead simply 
assume f = 0(1). Because our main focus in this paper concerns 
moderately to highly crowded systems, 0 = 0(1), for simplicity 
we set / = 1 throughout. Accordingly we work simply with 

* = ■£? (22) 


In the regime f <1, particles change their orientation many 
times during a typical interval between collisions. Therefore the 
dynamics during an interval 0 < t < t c between any two collisions 
will comprise only a proportionately much shorter regime of bal¬ 
listic motion at speed v a for times 0 < t < z 0 <C t c , followed by a 
longer diffusive regime for times t 0 < t < t c . In this way the “mi¬ 
croscopic” (pre-collision) dynamics that feed forward to inform 
the collisional dynamics are predominantly diffusive, with an ef¬ 
fective diffusion coefficient v\/D v . It might then be reasonable to 
expect that the phase diagram of the ABPs will mirror that of the 
PBPs with the scaled diffusion coefficient D = Dy/e of the passive 
system replaced its scaled counterpart v\y/D v £ for the active one. 

In contrast, in the regime £ » 1 the dynamics of the ABPs will 
be dominated by ballistic swimming over the entire time between 
collisions, with the slow angular reorientation being relatively 
unimportant in comparison. In this case it is less clear, upfront, 
what correspondence might be expected between the phase dia¬ 
gram of the ABPs with its equilibrium counterpart for the PBPs. 
As elaborated further below, however, we will indeed find some 


correspondence in the regime of low activity. 

This existence of two different regimes of the dimensionless 
inverse angular diffusion coefficient £ points to two different pos¬ 
sible choices for adimensionalising the swim speed v a . The first, 
which we expect to be the natural choice in the regime £ 1, 

takes into account both ballistic swimming and angular reorien¬ 
tation by recognising that v\/D r has the dimensions of a trans¬ 
lational diffusion coefficient then writing v a = v a y/D r £, by direct 
analogy with D = Dy/e for the PBPs. 

The second choice, which we expect to be the natural one in 
the regime £ ^> 1, considers only the ballistic component of the 
swimming dynamics, taking the ratio of v a with the characteristic 
speed e/Ry that a particle would exhibit when subject to a force 
of magnitude e/R in a solvent of drag coefficient y. This gives an 
adimensionalised swim speed v a = v^Ry/e that ignores the angu¬ 
larly diffusive component of the dynamics as characterised by A- 
It is nonetheless still seen as a possible counterpart of the D of 
PBPs by writing D = Dy/e = (D/R)Ry/e, and noting that D/R has 
dimensions of speed. We note that for £ = 1, v a = v a and the two 
methods of adimensionalising v a coincide, as required. 

With this preamble in mind, we now present results for the 
phase behaviour of ABPs in the plane of (v a ,0) for two values 
of the adimensional inverse diffusion coefficient £, one in each 
regime just described, with the swim speed v a in each case first 
adimensionalised in a manner suited to the value of £ according 
to the arguments just outlined. Where useful, we will also further 
present results in the plane of (v a , £) for different values of 0. 

We start in Fig. [6] by showing long-time snapshots of the sys¬ 
tem’s configuration on a grid of values of (v a , 0) for a low value of 
the scaled inverse diffusion coefficient £ = 0.2. Immediately ap¬ 
parent is that this non-equilibrium phase diagram displays strik¬ 
ing similarities to its equilibrium counterpart for the PBPs in the 
plane of (A0): recall Fig. [3] In particular, a G-S coexistence re¬ 
gion is seen for values of v a less then about 0.4. Above this a G-L 
coexistence region extends up towards what resembles a critical 
point somewhere in the region of 0 = 0.4 and v a = 0.7, which is 
quite close to that identified above for the PBP. For high area frac¬ 
tion, solid-like ordering is also again evident. For this small value 
of £, then, the main features of the equilibrium phase behaviour 
of PBPs transcribe to the non-equilibrium ABPs, with the scaled 
temperature D = k^T/e replaced by the scaled swim speed v a . 

Again as for the PBP, in a lower-left rectangle of the plane 
(v a , 0) that roughly corresponds to 0 < 0.25 and v a < 0.4, the sys¬ 
tem doesn’t reach steady state in any feasible run-time but instead 
shows a characteristic power-law coarsening process in which the 
typical domain size grows in time as h ~ . The data are well fit 

by the exponent /3 =2/3, as in the PBP above (data not shown). 

We have established, then, that for small £ the phase behaviour 
(and coarsening kinetics) of ABP in the plane of (v a , 0) closely re¬ 
sembles that of PBP in the plane of (A0). This correspondence 
is with hindsight perhaps unsurprising, given that the intercolli- 
sional dynamics are mainly diffusive in this regime, as in a passive 
system. 

We now turn to the opposite regime of large £, where v a is 
the most obvious choice for adimensionalising v a . Fig. [7] shows 
long-time snapshots of the system’s configuration on a grid of val- 
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Fig. 6 Active Brownian particles, £ = 0.2. Snapshots of the system’s configuration on a grid of values of the scaled active swim speed v a and effective 
area fraction 0. Each snapshot is taken at a long time 3 x \0 3 v a /R after the system was initialised at t = 0 in a random state as described in Sec. [ 4 ] 
Note that the scale is nonlinear at the largest v a and at low 0. 
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Fig. 7 Active Brownian particles, f = 200.0. Snapshots of the system’s configuration on a grid of values of the scaled active swim speed v a and 
effective area fraction 0. Each snapshot is taken at a long time 3 x 10 3 v a /R after the system was initialised at t = 0 in a random state as described in 
Sec.[4] Note that the scale is nonlinear at the largest v a and at low 
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Fig. 8 Active Brownian Particles. Area fraction fluctuations 8$ for a 
small (upper) and large (lower) value of £. Quasi-equilibrium phase 
separation is seen at low scaled active swim speed v a for both values of 
£. Strongly non-equilibrium MIPS is seen at high v a only for large £. 
Number of boxes G = 5. 


ues of (v a ,0) for £ = 200.0. The situation here is clearly more 
complicated than for the low value of £, with the phase diagram 
differing more markedly from the equilibrium case. This is to be 
expected, consistent with the fact that underlying “microscopic” 
dynamics that feeds forward to inform the collisional dynamics is 
predominantly ballistic, rather than diffusive as it was in the case 
either of PBPs or ABPs for small £. 

Nonetheless, at low swim speeds the phase behaviour does ap¬ 
pear rather similar to that of PBPs, and so also of ABPs at small 
values of £. Accordingly we suggest that this again corresponds 
to a regime of phase separation analogous to that of the equilib¬ 
rium PBPs, arising due to the dominance of the Lennard-Jones 
attraction over the active swimming for small v a . However we 
repeat that the physics is subtly different here: unlike the equilib¬ 
rium case the intercollisional dynamics is predominantly ballistic 
rather than diffusive. In view of this, the fact that we do still see 
a regime of attraction-dominated quasi-equilibrium phase separa¬ 
tion at low v a is a non-trivial result. On this basis, one might tenta¬ 
tively claim that the ABPs’ scaled swim speed has some properties 
resembling those of a non-equilibrium effective temperature. 

At high activity (large v a ), however, the phase behaviour for 
this large value of £ is markedly different from the equilibrium 
case. In particular, a second regime of phase separation is evident 


Fig. 9 Active Brownian particles: area fraction fluctuations 8$ in the 
8 - £ plane. The black line corresponds to the curve v a £/(l + £) = 1. 
See Eqn.[23lin the main text. Number of boxes G = 5. 


for values of 0 exceeding about 0.4. Clearly, this could not have 
been predicted on the grounds of any direct analogy with the 
equilibrium phase diagram. How then can it be understood? 

For these large values of v a we expect the strength of swim¬ 
ming to easily overcome any attraction associated with the LJ 
potential, such that active LJ disks behave effectively as active 
hard disks in this limit v a —>• «>. With this in mind, we recall pre¬ 
vious work for hard disk ABPs that demonstrated a phenomenon 
known as “motility induced phase separation” (MIPS) 21 JHI This 
arises via a positive feedback mechanism in which particles (a) 
accumulate in regions where they move more slowly, then (b) 
further slow down in regions of accumulation, being impeded by 
overcrowding from other particles. Note that although part (a) of 
this mechanism may seem intuitively obvious - think of window 
shoppers accumulating where they linger near a more interesting 
display - it is actually a highly non-equilibrium phenomenon: in 
equilibrium colloids the local density depends only on the static 
potential and not on any dynamical rules. Informed by this, we 
interpret this second region of phase separation, seen at high v a 
for large £, as one of MIPS. 

The absence of MIPS at high v a for small values of £, as in 
Fig. [6] is then further understoo d! 28 ! 61 ! by realising that part (b) 
of the feedback mechanism just described contains a mean field 
assumption that fails for small £: in order to experience a slowing 
down of its run speed, a particle must collide many times during 
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any such run before changing direction. This assumption breaks 
down for small values £, where particles reorient very quickly on 
the timescales of interparticle collisions. 

This difference between the phase behaviour of ABPs for small 
and large values of £ is characterised at a glance by colour maps 
of the particle number fluctuations in Fig. [8] MIPS is apparent for 
large £ (bottom panel) for 0 ~ 0.5 and v a ~ 1.0; and absent for 
small £ (top panel). 

The dependence on £ is shown more fully in Fig. [9] For the 
small area fraction 0 = 0.125, no MIPS is expected for any value 
of £ (recall Fig. [8]) and we see only quasi-equilibrium phase sep¬ 
aration driven by the Lennard-Jones attractions. As can be seen, 
this arises for v a < 0(1) at large £, consistent with the expectation 
that v a = v a Ry/£ is the relevant scaling variable in this regime. In 
contrast, for low £ it arises instead for v a < 0(1), consistent with 
the expectation that v a = v\y/D Y £ is the relevant scaling variable 
in this regime. Noting that v a = v a £, we find finally that 


V~a£ 

1 + C 


< 1 


(23) 


gives a good fit to the regime of attraction-dominated phase coex¬ 
istence for all values of £, recovering v a < 1 for £ <c 1 and v a < 1 
for £ > 1. 

For the larger area fraction 0 = 0.75 in Fig.[9]we find the same 
picture with regards attraction-dominated phase separation, with 
the transition again given by Eqn. [23] In addition to this we find 
a regime of MIPS for large values of £ ~ 10.0 and v a ~ 1.0. Inter¬ 
estingly, these regimes of equilibrium-like attraction-dominated 
phase separation and MIPS appear to merge with each other 
smoothly, with a gradual cross over regime in which both mecha¬ 
nisms contribute. We note that re-entrance between two regimes 
of phase separation was also observed in simulations of attrac¬ 
tive ABPs with the conventional LJ potential in Ref. 47 . However 
that work considered only the regime of small rotational diffusion 
coefficient (of thermal origin, and so related to the Brownian dif¬ 
fusion coefficient via a factor 3/ R 2 ) and therefore necessarily saw 
a regime with MIPS. The present study goes beyond the findings 
of 47 in showing MIPS to be suppressed with increasing rotational 
diffusion. (RBll/15)In Ref.SSj, a crossover from equilibrium-like 
phase behaviour to a non-equilibrium percolating network phase 
was seen in a system of active Brownian Lennard Jones particles 
as function of decreasing rotational diffusion coefficient. The per¬ 
colating network phase in that study might well be the counter¬ 
part of the MIPS seen here, though further work would be needed 
to elucidate further any connection. 


6.3 Hydrodynamic squirmers 

In the previous section, we explored the phase behaviour of active 
Brownian particles (ABPs), which combine ballistic swimming 
with stochastic angular reorientation that is intended to mimic, 
in continuous-time counterpart, discrete tumbling events in some 
species of bacteria, or scattering off other particles due to hydro- 
dynamic interactions. Hydrodynamic interactions are not, how¬ 
ever, explicitly accounted for in ABPs. In this section we turn to 
address that shortcoming by studying a suspension of squirmers 


that do interact hydrodynamically. 

We start with a recap of the parameters that characterise 
squirmers and the appropriate adimensional form of these, as first 
introduced in Sec. [3] above. Indeed to set these in a proper con¬ 
text, we recall them alongside the (by now) familiar parameters 
of the PBP and ABPs. 

As always we take as fixed the number of particles N = 242 
and the lengthscales <Jh and Olj that set the size of the poten¬ 
tial shell relative to the hydrodynamic radius R. For the PBPs the 
remaining five parameters are then the particle radius R, the po¬ 
tential energy e, the drag coefficient y, the area fraction 0 and 
the thermal translational diffusion coefficient D. Choosing units 
of mass, length and time left two remaining dimensionless pa¬ 
rameters: D = Dy/e = k^T/s and 0. In Sec. 6.1 we explored the 
phase behaviour of our modified Lennard-Jones particles in this 
plane of (0,0). For the ABPs the remaining parameters were the 
particle radius R, the potential energy e, the drag coefficient y, 
the area fraction 0, the swim speed v a and the rotational diffu¬ 
sion coefficient D v . In dimensionless form these gave three pa¬ 
rameters, which we explored numerically in the previous section. 
Two of these were the area fraction 0 and the scaled inverse ro¬ 
tational diffusion coefficient £ = v a /0 r /?, which characterises the 
time taken for the angular decorrelation of swim direction rela¬ 
tive to the typical time interval between particle collisions. The 
third was the swim speed v a , for which we argued v a = v\y/D v £ 
to be an appropriate adimensional form in the regime of small 
£ where the inter-collisional dynamics are diffusion dominated, 
and v a = v^Ry/e at high £, where the inter-collisional dynamics 
are predominantly ballistic. Given that v a = v a £ these two coin¬ 
cide for values of £ = 0(1), as desired. 

For the hydrodynamic squirmers to which we now turn the rel¬ 
evant parameters are the particle radius R, the potential depth 
£, the solvent viscosity q, the area fraction 0, the swim speed 
v a = B i/2 and the relative stresslet strength /3 = B 2 /B 1 . In dimen¬ 
sionless form we then have the scaled swim speed v a = v a q R 2 /£, 
the area fraction 0, and the stresslet strength j 3. We shall show in 
what follows that the effects of j8 on phase behaviour are in fact 
relatively mild. Accordingly we shall focus first, and mainly, on 
phase behaviour as a function of v a and 0, before returning at the 
end of the section to comment on the effects of varying /3. 

Excluding j8 for the moment, then, we note the number of im¬ 
portant parameters (v a , 0) for the squirmers to be one fewer than 
those for the ABPs, (v a , 0, £) (or (v a , 0, £)). The reason for this is 
as follows. Whereas the ABPs have externally imposed reorienta¬ 
tion dynamics with a rotational diffusion coefficient D r , adimen- 
sionalised as £, the squirmers have no reorientation dynamics im¬ 
posed upfront. Instead, reorientation events emerge naturally as 
a result of the hydrodynamic interactions that determine the way 
in which the squirmers scatter off each other (at low area frac¬ 
tion) or slither round each other (at high area fraction). In con¬ 
sequence we expect an effective ratio £ e ff = t 0 /t c of reorientation 
time to scattering time to emerge naturally from the squirmers’ 
many body hydrodynamics, instead of existing as a free parame¬ 
ter to be imposed at the outset in the simulations. 

Having identified the relevant parameters, a natural question is 
then whether the phase behaviour of the squirmers in the plane of 


This journal is ©The Royal Society of Chemistry [year] 


Journal Name, [year], [vol.],i -[ 23 ] | 15 




v a 

0.275 

0.25 

0.225 

0.2 

0.175 

0.15 

0.125 

0.1 

0.075 

0.05 

0.025 



0.016 0.031 0.063 0.094 0.125 0.25 0.375 0.50 0.625 0.75 0.875 <p 


Fig. 10 Hydrodynamic squirmers, j8 = 0.0. Snapshots of the system’s configuration on a grid of values of the scaled active swim speed v a and 
effective area fraction </>. Each snapshot is taken at a long time after the system was initialised at t = 0 in a random state as described in Sec.[4] In the 
rectangle v a < 0.15 and 0 < 0.1 the system fails to reach steady state in any feasible runtime and instead shows slow domain coarsening. All other 
snapshots shown are from a system in statistically steady state. Note that the scale is nonlinear at the largest v a and at low </). 
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(v a = v a T 7 /? 2 /£, 0 ), with whatever value of £ e ff emerges naturally 
from the hydrodynamic simulations, mimics that of the ABPs in 
the plane (v a = v a y/?/e,0) for a corresponding value of £, as im¬ 
posed for the ABPs. One can of course also further ask whether it 
further mimics that of PBPs in the plane of (Z), 0). 

Clearly, the dimensionless forms for v a identified above sug¬ 
gest that we might look for a possible correspondence between 
ABPs and squirmers (after matching the imposed £ for the ABPs 
with the emergent Ceff for the squirmers) upon setting the ABPs’ 
drag coefficient y= rjR, where q is the solvent viscosity for the 
squirmer simulations. However it should be noted that this ex¬ 
pression should also contain a prefactor, which is unknown. The 
specific case of a sphere dragged by a force monopole through an 
infinite domain of fluid in 3D gives y = 6nriR. However we are 
concerned here neither with force monopoles nor with 3D and 
the prefactor remains unknown, even if we might reasonably ex¬ 
pect it to be 0(1). With this preamble in mind we now present 
our results. We start by taking a fixed /3 = 0.0 and shall return 
at the end of the section to comment on the effect of varying j3. 
Fig. [To] shows snapshots of the system’s configuration on a grid 
of values of (v a , 0). Each snapshot corresponds to the longest run 
time that is feasibly accessible computationally. Most represent a 
statistically steady state, except in a region towards the lower left 
of the (v a ,0) plane, typically for values of v a < 0.15 and 0 <0.1. 
Here the coarsening dynamics of the aggregated domains, which 
we shall discuss in more detail below, is sufficiently slow that the 
system never attains steady state in this regime in any reasonable 
run time. Despite this we are confident that running for even 
longer times would not affect the qualitative features of the sys¬ 
tem’s state, apart from giving larger domains. 

Some quantitative measures of this phase diagram are shown 
in Fig. EH Panel a) gives the mean number of particles per 
cluster h, normalised by the total number of particles N to give 
h* = h/N. Recall that a single system-spanning condensate would 
have h* = 1, while a gas phase with no clusters would have 
n* = 1 /N, which is small for large N. Panel b) shows the particle 
number fluctuations 8^/VR measured by dividing the simulation 
box into G x G subboxes where G = 5. Subpanel c) shows the 
hexatic order parameter 'Pg. Recall that regions of high h* and 
8n/Vn indicate clustering or phase separation, while regions of 
high 'Pg indicate a high degree of solid-like crystalline ordering. 

Comparing these snapshots and colourmaps for the squirmers 
(Figs. [lO] and [TT] > with their counterparts for the PBPs (Figs. [3] 
an d|4]l reveals that the overall phase behaviour of the squirmers 
loosely mirrors that of the PBPs, with the scaled swim speed v a 
playing a role somewhat analogous to that of the scaled diffu¬ 
sion coefficient D of the PBPs. In particular, attraction-dominated 
phase separation and hexatic ordering appears evident at low v a 
(corresponding to low D for the PBPs), giving way to disorder and 
homogeneity at high v a (high D for the PBPs) where the particles’ 
activity (thermal jostling for the PBPs) is sufficient to overcome 
the attractive effects of the potential. Solid-like ordering is also 
evident in a column of values of v a for densely crowded systems 
with 0 > 0.75. 

However there are also important differences. For example, in 
the squirmers hexatic ordering (Fig. [TT] :) is only strongly appar- 
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Fig. 11 Phase diagram of hydrodynamic squirmers mapped using (a) 
normalised mean number of particles per cluster n* = n/N, (b) particle 
number fluctuations S N /Vti for G = 5 and c) hexatic order parameter TV 


ent at the lowest values of v a for which appreciable area fraction 
fluctuations are present (Fig. [TTJ)), reflecting the fact that the ev¬ 
idently phase separated snapshots in Fig. [lO] only show solid-like 
ordering at the lowest v a . In contrast, for the PBPs hexatic order 
(Fig. [ 4 ]:) persists to higher values of D relative to those for which 
appreciable phase separation is evident (Fig. ER The same is 
true of the column of hexatic ordering seen at high 0: this grad¬ 
ually fades with increasing v a in the squirmers, whereas it appar¬ 
ently persists to indefinitely high D for the PBPs in Fig. &• 

This gradual decay of hexatic order with increasing v a suggests 
that the attraction-dominated phase ordering is disrupted more 
readily by the ballistic intercollisional swimming of the squirm- 
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ers with increasing v a than it is by the diffusive intercollisional 
dynamics of the PBP with increasing D. This is consistent with 
the fact that these coherently moving “living crystals” seen in the 
regime of clustering at low v a in Fig. [To] are much more dynamic 
and mobile entities than their equilibrium counterparts in Fig. [3] 
The mean squared displacement (MSD) of particles in these clus¬ 
ters (data not shown) have an early-time ballistic regime crossing 
over into a later-time diffusive regime, even in crowded systems. 
In contrast, the counterpart MSDs for crowded PBPs at low D have 
a slower-than-diffusive (“caged”) regime at early times crossing 
over into diffusive dynamics later on (data not shown). Movies 
of these clusters can be found online^ and are seen to resemble 
those observed experimentally in active colloids in Refs.ESHIESl 
named as “living crystals”. 

Another notable difference is that stringlike clusters remain 
clearly evident in the squirmer snapshots even at high v a (Fig.fTo]), 
compared with the more complete return to homogeneity at high 
D in the PBPs (Fig. [ 3 }. This is consistent with the existence of a 
larger window between the lines of h* = 1 /50 and 1/10 in Fig.[TTfc 
compared to its equilibrium counterpart in Fig. Eh indicating a 
larger regime of clustering in the squirmers than in the PBPs. 
Likewise the area fraction fluctuations decay more slowly as a 
function of v a in Fig. Y\_ than they do as a function of D in Fig. [ 4 ] 

Clearly, this slight clustering of squirmers in the regime of high 
v a has no counterpart in the equilibrium phase diagram of the 
PBPs at high D. It is, however, reminiscent of the phase behaviour 
of ABPs for low to moderate values of £, where a slight clustering 
arises at high v a , which we interpreted as a precursor to the onset 
of true MIPS at high £. 

Tentatively, then, we interpret this tendency of the squirm¬ 
ers to form stringlike clusters at high v a as a counterpart of the 
slight clustering seen in the ABPs for modest imposed values of 
£ = 0(1). Recalling that £ corresponds to the ratio of reorien¬ 
tation time and (estimated) collision time t 0 /t c that is imposed 
upfront in the ABP simulations, if this interpretation is correct we 
might then reasonably anticipate that the corresponding effec¬ 
tive ratio £ e ff = t 0 /t c that emerges from the squirmer simulations 
should likewise be 0(1). This prediction of an emergent £ = 0(1) 
was indeed confirmed (for the case of hard disk squirmers) in 
Ref.ES. it is furthermore consistent with earlier observations that 
whenever two squirmers scatter off each other hydrodynamically 
they also suffer an 0(1) change in their swim directions SO. To re¬ 
inforce this argument we extracted from our squirmer simulations 
this emergent effective £ e ff, measuring the characteristic times t 0 
and t c in the way defined in Ref. ISO. As can be seen in Fig. 


12 


we 

indeed recover 0(1) (or smaller) values for the effective emergent 
£ e ff for all the mid-high area fractions, which would be the ones 
potentially inside the MIPS regime. 

We therefore suggest that a slight tendency to clustering seen at 
high v a in the squirmers is the signature of a nearby MIPS, the full 
effects of which have been mitigated by the fact that the effective 
emergent value of £ e ff for the squirmers is only 0(1), rather than 
the larger value of £ needed to see full MIPS in the ABPs. The 
underlying physics here is that, as just noted, squirmers change 
their orientations each time they scatter off each other hydrody¬ 
namically, thereby giving an effective £ e ff = 0(1). In contrast for 


the ABPs we can artificially tune £ to be arbitrarily large simply 
by setting a smaller D v in the simulations, thereby inducing MIPS. 
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Fig. 12 Ratio between the characteristic reorientation time t d and 
collision time t c for polar squirmers, giving the emergent effective ratio 

Ceff- 


As noted above, this suppression of MIPS by hydrodynamics 
was first discussed in the context of purely repulsive particles in 
Ref.ES. a contribution of the present work has been to estab¬ 
lish that MIPS is also suppressed by hydrodynamics in systems 
of attractive particles, despite any intuition one might have held 
upfront that attractive interactions might enhance the crowding 
mechanism of MIPS. The situation could however be different in 
systems such as in Janus particles!® where the interaction poten¬ 
tial is orientation-dependent and might be expected to increase 
the angular decorrelation time t g , thereby increasing the ten¬ 
dency towards MIPS. 

Although we have demonstrated the suppression of MIPS by 
hydrodynamics only in this athermal limit of D v = 0 and D = 0 for 
which we simulate the squirmers, it was in this limit that MIPS 
was found to be most pronounced in systems of ABPs without hy¬ 
drodynamics in Ref. 27 . This gives a strong indication that MIPS 
should remain suppressed even at non-zero temperature, D r ^ 0 
and D / 0. Indeed, any extra thermal contribution to D v would be 
expected only further to decrease the ratio of £ e ff = t 0 /t c , thereby 
further suppressing MIPS. The same argument would apply to 
the explicit inclusion of (athermal) tumbling events in the simu¬ 
lations. 

To summarise, we suggest that the phase behaviour of attrac¬ 
tive squirmers in the plane of (v a ,0) can be understood in terms 
of two different mechanisms that smoothly blend into each other. 
The first is a regime of attraction-dominated phase separation 
at small v a , which we suggest is analogous to the equilibrium 
phase separation seen in PBPs at small D, although more easily 
disturbed by the ballistic intercollisional dynamics of the squirm¬ 
ers as v a increases. The second is a regime of string-like particle 
clustering at high v a , which we have interpreted in terms of a 
nearby MIPS that has been mitigated by the fact that £ e ff = t 0 /t c 
is automatically 0(1) for squirmers, which reorient each time they 
scatter off each other. 

So far, we have presented results for the hydrodynamic squirm¬ 
ers all for a single value of the stresslet parameter /3 = B 2 /B 1 = 
0.0. We shall now briefly consider the effects on phase behaviour 
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Fig. 13 Hydrodynamic squirmers. Area fraction fluctuations <50 
obtained on dividing the simulation box into GxG subboxes with G = 5. 

of varying /3. Recall that a single particle in an infinite medium 
swims at a speed v a = B\/2, which is set only by B\. The param¬ 
eter #2 instead measures an additional apolar contribution to the 
flow field that establishes round a (single, isolated) swimmer, and 
gives no contribution to the single-particle swim speed. In terms 
of the ratio j8 = B 2 /B i, therefore, the regime |/3| <C 1 corresponds 
to strongly polar swimmers while |j8| 1 corresponds to apolar 

“shakers”. 

Although B 2 has no effect on the swim speed of a single particle 
in an infinite medium, in the case of swimmers crowded together 
we might anticipate that the additional apolar flows governed by 
B 2 do lead to higher swim speeds for all the squirmers collectively, 



1 10 100 


tv a /R 

Fig. 14 Hydrodynamic squirmers, p = 0 .0. Temporal evolution of the 
number of particles per cluster, n*(t) (minus the initial contribution 
n*(t = 0)). Time is shown in units of the characteristic time v a /R that an 
isolated particle would take to displace one particle radius. The dotted 
line shows the power law r 2 / 3 . 

as the result of hydrodynamic interactions between them: the 
apolar contribution to the flow field round one squirmer disturbs 
neighbouring squirmers, causing them to move slightly faster. 

In view of this, we suggest that the adimensionalised measure 
of activity v a on the vertical axis of the phase diagram should now 
be generalised from the bare quantity v a = B\r\R 2 /2e = v&rjR 2 /e 
used so far to a quantity that more fairly reflects the additional 
contribution of B 2 to the overall collective activity. Accordingly 
we combine B\ and B 2 as \Jb 2 + a 2 B 2 , and define an updated 
adimensional v a = v a ^/l + a 2 fi 2 r\R 2 /e. (This reduces as required 
to the quantity v a = v a q R 2 /e used to represent the /3 = 0 results 
above.) A key question then is whether the squirmers’ phase be¬ 
haviour in the plane of (this newly defined) v a and 0 depends in 
any significant way on /3. 

The parameter a, which is defined by the formulae just given, 
deserves some discussion. In general it should depend on the area 
fraction 0, because the degree to which the apolar components of 
the flow field round each particle will cause neighbouring parti¬ 
cles to move via hydrodynamic interactions will depend also on 
the particle separation and so on 0, consistent with the fact that 
B 2 has no effect on swimming dynamics in the single particle limit 
0 —0, where only B\ matters. However we shall ignore this com¬ 
plication in what follows and simply set a constant throughout. 

Fig. [13] shows the phase behaviour of the squirmers, as mea¬ 
sured by the particle number fluctuations, in this plane of v a = 
v a \/l + oc 2 p 2 riR 2 /e and 0 for a value a = 0.1. For this particu¬ 
lar chosen value of a, the main features of the phase diagram 
are preserved in both magnitude and location even as /3 is varied 
between 0,5,°°. Although all the panels shown in this figure are 
for pullers (/3 >0), very similar results (not shown) are found for 
pushers (/3 < 0). Provided the data is represented in this way, we 
conclude that the stresslet parameter /3 is relatively unimportant 
in determining the clustering and phase behaviour of active par¬ 
ticles with hydrodynamics, and that disklike “pushers”, “pullers” 
and “shakers” should all show similar behaviour. The small value 
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of a needed to ensure this - for other values the region of strong 
fluctuations shifts in vertical position with changing /3 - suggests 
that #2 has a relatively large effect compared with B\ in deter¬ 
mining the overall magnitude of effective activity. This is consis¬ 
tent with the notion that B\ acts even in the limit of zero volume 
fraction whereas B 2 affects the activity levels through interactions 
only. Accordingly the effect of B 2 should scale as the inverse parti¬ 
cle separation, which is large for the relatively crowded, clustered 
regimes of interest here. 

Finally, we turn from the properties of the statistically steady 
state that forms a long time after the system was randomly ini¬ 
tialised at time t = 0, to comment briefly on the dynamics of the 
aggregated domains that slowly coarsen at low v a . See Fig. [14] 
which shows that the number of particles per domain h*{t) ap¬ 
pears to coarsen as a power law t 2 / 3 , corresponding to a char¬ 
acteristic linear domain radius t 1 / 3 , which we note is the same 
power law as for the coarsening dynamics of an undriven system 
in a diffusion dominated regime^. Given that the squirmers have 
hydrodynamic interactions and ballistic rather than diffusive in- 
tercollisional dynamics, this correspondence is a nontrivial and 
perhaps even surprising result. 

7 Conclusions 

Motivated by recent experiments reporting clustering in active 
suspensions^® 32 , we have studied with full hydrodynamics the 
phase behaviour of active particles with an interaction potential 
that has an attractive Lennard-Jones (LJ) component, together 
with a steep repulsive contribution that diverges at particle con¬ 
tact r 2R+, where R is the particle radius. 

As a preliminary step we mapped out the equilibrium phase di¬ 
agram of passive particles subject to the same interaction poten¬ 
tial, finding phase behaviour broadly as expected for a LJ system, 
with regions of G-S, G-L and L-S phase coexistence at low tem¬ 
peratures. However compared to the LJ potential used more con¬ 
ventionally in the literature, the repulsive component of which di¬ 
verges only as r —»0, we found crystalline order to set in a slightly 
lower packing fractions, consistent with the greater excluded vol¬ 
ume effects from the additional steeply repulsive contribution to 
our modified potential as r —>> 2 R. This preliminary study also es¬ 
tablished that the main features of phase behaviour can indeed be 
mapped out with only a small number of particles, N = 242, which 
is the maximum feasible system size in the much more costly hy¬ 
drodynamic simulations. 

As a second preliminary step, we then studied active Brownian 
particles (ABPs) subject to the same potential, but without hydro¬ 
dynamics. As in previous studies of ABPs with the conventional 
Lennard-Jones potential 47 !, we observed two distinct regimes 
of phase separation. For low activities a regime of attraction- 
induced phase separation is seen, strikingly ressembling that at 
low temperatures in the equilibrium phase diagram, with a return 
to more homegeneous states for higher activities. In this sense, 
activity acts in a manner somewhat analogous to the temperature 
of a passive system. For ABPs with a small rotational diffusion 
coefficient, however, we found at high activities a second regime 
of phase separation that has no counterpart in the passive system 
and is instead a purely non-equilibrium phenomenon known as 


motility induced phase separation (MIPS) 2HH1. 

Motivated by the important effects of hydrodynamics in many 
active systems, we proceeded finally to simulate a suspension of 
hydrodynamic squirmers subject to the same (modified) LJ poten¬ 
tial. The central results of this hydrodynamic study are fourfold. 

First, we showed that activity again acts in a manner somewhat 
analogous to an effective temperature, in the sense that a regime 
of equilibrium-like attraction-induced phase separation is seen at 
low activities (with the activity suitably adimensionalised by the 
strength of the attractive potential), with a return to much more 
homogenised states at high activity. 

Second, at high values of the (scaled) activity we find only 
weak stringlike clustering, with no evidence for motility induced 
phase separation (MIPS). By careful comparison with the ABPs, 
we interpreted this weak clustering as the signature of a nearby 
MIPS that has been largely suppressed by hydrodynamic interac¬ 
tions according to the arguments in RefGS. 

In view of this suppression, we tentatively suggest that cluster¬ 
ing phenomena seen experimentally in active colloids in which 
hydrodynamics is important are more likely to be in the regime of 
equilibrium-like attraction-dominated phase separation than the 
regime of purely non-equilibrium MIPS. Indeed, as noted above, 
the presence of a small phoretic attraction can seldom be ruled 
out in synthetic colloids. It is also interesting to recall that clus¬ 
tering is seldom seen in bacterial suspensions 33 “S3 at least in the 
absence of an obvious depletion interactionISEEl 

Third, we showed that the regime of phase separation seen at 
low activity, despite being broadly analogous to its equilibrium 
counterpart in the passive systems, also shows important differ¬ 
ences. In particular we find a lower than expected degree of hex- 
atic ordering, compared to the degree of clustering, consistent 
with the more dynamic nature of so these “living crystals” B°^ 32 i 
Activity is therefore seen to disrupt the level of translational or¬ 
dering compared to that in the corresponding equilibrium phase 
diagram. 

Fourth, we showed that an extended definition of the rescaled 
activity parameter allows us to recover strikingly similar phase di¬ 
agrams across the full range of the stresslet parameter j8, suggest¬ 
ing that “pushers”, “pullers” and “shakers” all behave in a similar 
manner with regards their aggregation and phase behaviour. 

In future work it would be interesting to consider more fully 
the consequences of dimensionality on these effects. Recall that 
the system simulated in this work concerns a 2D layer of ac¬ 
tive particles with 2D hydrodynamic interactions between them. 
While we are confident, on the basis of checks reported earlier 
in Ref.^U, that the same effects would also obtain for a 2D layer 
of particles with 3D hydrodynamic interactions, it would clearly 
be interesting to consider the case of 3D packings with 3D hydro¬ 
dynamics. Indeed, loose aggregates moving coherently for long 
times were observed previously for a 3D suspension of purely re¬ 
pulsive squirmers with 3D hydrodynamics in Refill. In Ref.^H 
a strongly confined 2D monolayer of purely repulsive squirm¬ 
ers was observed to show strong phase separation. However it 
should be noted that an important effect of strong confinement 
will be to screen hydrodynamic interactions. In this sense, there¬ 
fore, the observation of strong phase separation in Ref. 61 is per- 
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haps not surprising. Far more puzzlingly, phase separation was 
apparently observed in an unconfined 2D monolayer of purely re¬ 
pulsive squirmers with fully 3D hydrodynamics in RefJSl though 
at an area fraction 0 = 0.1 that seems far too low for MIPS to be 
implicated. 

Asymmetric particles could be considered in order to mimic 
chemotactic swimmers, with their rich phase diagramsEIKO A 
closer connection with experiments on phoretic active colloids 
might be established by solving the dynamics of a chemical con¬ 
centration field in the fluid between the particles, with a suitable 
boundary condition at the particle surfaces, to account for diffu- 
siophoresis. 

Finally, in view of the important role played by the ratio £ of 
the angular reorientation time to the interparticle collision time, 
it would be interesting to simulate elongated particles with full 
hydrodynamics. One might anticipate that particle elongation 
would significantly affect the angular reorientation time, perhaps 
leading to an effective emergent value of £, even with hydrody¬ 
namics, that would predispose a return to MIPS at high values of 
activity. It remains an open challenge to determine the link, if any, 
between that possible mechanism for MIPS in elongated particles 
and the one originally put forward in Refs .USES i n the context of 
active liquids crystals. 
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